Revisiting the Phase Transition of Spin-1/2 Heisenberg Model with a Spatially 
Staggered Anisotropy on the Square Lattice 
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Puzzled by the indication of a new critical theory for the spin-1/2 Heisenberg model with a 
spatially staggered anisotropy on the square lattice as suggested in [J, we re-investigate the phase 
transition of this model induced by dimerization using first principle Monte Carlo simulations. 
We focus on studying the finite-size scaling of p s iL and p&L, where L stands for the spatial box 
size used in the simulations and p s i with i £ {1,2} is the spin-stiffness in i-direction. From our 
Monte Carlo data, we find that p S 2L suffers a much less severe correction compared to that of p 3 \L. 
Fherefore p S 2L is a better quantity than p s \L for finite-size scaling analysis concerning the limitation 
of the availability of large volumes data in our study. Further, motivated by the so-called cubical 
regime in magnon chiral perturbation theory, we additionally perform a finite-size scaling analysis 
on our Monte Carlo data with the assumption that the ratio of spatial winding numbers squared 
is fixed through all simulations. As a result, the physical shape of the system remains fixed in our 
calculations. The validity of this new idea is confirmed by studying the phase transition driven by 
spatial anisotropy for the ladder anisotropic Heisenberg model. With this new strategy, even from 
PsiL which receives the most serious correction among the observables considered in this study, we 
arrive at a value for the critical exponent v which is consistent with the expected 0(3) value by 
using only up to L — 64 data points. 

PACS numbers: 



I. INTRODUCTION 



Heiscnberg-type models have been studied in great de- 
tail during the last twenty years because of their phe- 
nomenological importance. For example, it is believed 
that the spin-1/2 Heisenberg model on the square lattice 
is the correct model for understanding the undoped pre- 
cursors of high T c cuprates (undoped antiferromagnets) . 
Further, due to the availability of efficient Monte Carlo 
algorithms as well as the increasing power of computing 
resources, properties of undoped antiferromagnets on ge- 
ometrically non-frustrated lattices have been determined 
to unprecedented accuracy 2 8]. For instance, using a 
loop algorithm, the low-energy parameters of the spin- 
1/2 Heisenberg model on the square lattice are calculated 
very precisely and are in quantitative agreement with the 
experimental results @ . Despite being well studied, sev- 
eral recent numerical investigation of anisotropic Heisen- 
berg models have led to unexpected results [l], [Io|, [H| ■ 
In particular, Monte Carlo evidence indicates that the 
anisotropic Heisenberg model with staggered arrange- 
ment of the antiferromagnetic couplings may belong to a 
new universality class, in contradiction to the theoretical 
0(3) universality prediction [l[. For example, while the 
most accurate Monte Carlo value for the critical exponent 
v in the 0(3) universality class is given by v = 0.7112(5) 
fl2j , the corresponding v determined in [l| is shown to be 
v = 0.689(5). Although subtlety of calculating the crit- 
ical exponent v from performing finite-size scaling anal- 
ysis is demonstrated for a similar anisotropic Heisenberg 
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model on the honeycomb lattice [13[ , the discrepancy be- 
tween v = 0.689(5) and v = 0.7112(5) observed in [H0J] 
remains to be understood. 

In order to clarify this issue further, we have simulated 
the spin-1/2 Heisenberg model with a spatially staggered 
anisotropy on the square lattice. Further, we choose to 
analyze the finite-size scaling of the observables p s \L and 
p S 2L, where L refers to the box size used in the sim- 
ulations and p s i with i £ {1,2} is the spin stiffness in 
i-direction. The reason for choosing p s \L and p s iL is 
twofold. First of all, these two observables can be cal- 
culated to a very high accuracy using loop algorithms. 
Secondly, one can measure p s \ and p S 2 separately. In 
practice, one would naturally use p s which is the aver- 
age of pgi and p S 2 for the data analysis. However for 
the model considered here, we find it is useful to analyze 
both the data of p s \ and p S 2 because studying p s i and 
p S 2 individually might reveal the impact of anisotropy 
on the system. Surprisingly, as we will show later, the 
observable p S 2L receives a much less severe correction 
than p s \L does. Hence p S 2L is a better observable than 
p s \L (or p s L) for finite-size scaling analysis concerning 
the limitation of the availability of large volumes data 
in this study. Further, motivated by the so-called cu- 
bical regime in magnon chiral perturbation theory, we 
have performed an additional finite-size scaling analy- 
sis on Ps\L with the assumption that the ratio of spa- 
tial winding numbers squared is fixed in all our Monte- 
Carlo simulations. In other word, we keep the physical 
shape of the system fixed in the additional analysis of 
finite-size scaling. The validity of this new idea is con- 
firmed by studying the phase transition driven by spatial 
anisotropy for the ladder anisotropic Heisenberg model, 
namely the critical point and critical exponent v for this 
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FIG. 1: 

study. 



The anisotropic Heisenberg model considered in this 



phase transition we obtain by fixing the ratio of spatial 
winding numbers squared are consistent with the known 
results in the literature. Remarkably, combining the idea 
of fixing the ratio of spatial winding numbers squared 
in the simulations and finite-size scaling analysis, unlike 
the unconventional value for v observed in [l[ , even from 
p s \L which suffers a very serious correction, we arrive at 
a value for v which is consistent with that of 0(3) by 
using only up to L = 64 data points. 

This paper is organized as follows. In section [III the 
anisotropic Heisenberg model and the relevant observ- 
ables studied in this work are briefly described. Section 
IIIII contains our numerical results. In particular, the cor- 
responding critical point as well as the critical exponent 
v are determined by fitting the numerical data to their 
predicted critical behavior near the transition. Finally, 
we conclude our study in section HVl 
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FIG. 2: Monte Carlo data of p s L (upper panel) and p S 2L 
(lower panel) as functions of the parameter J'/J. 



critical point as well as the critical exponent v with high 
precision. 



II. MICROSCOPIC MODEL AND 
CORRESPONDING OBSERVABLES 

The Heisenberg model considered in this study is de- 
fined by the Hamilton operator 

H = ]T JS X ■ S y + J' S*> ■ Sy, (1) 

(xy) {x'y') 

where J and J' are antiferromagnetic exchange couplings 
connecting nearest neighbor spins (xy) and (x'y'), re- 
spectively Figure 1 illustrates the Heisenberg model de- 
scribed by Eq. ([1]). To study the critical behavior of this 
anisotropic Heisenberg model near the transition driven 
by spatial anisotropy, in particular to determine the crit- 
ical point as well as the critical exponent u, the spin 
stiffnesses in 1- and 2-dircctions which arc defined by 

Psi = ^(Wl), (2) 

are measured in our simulations. Here j3 is the inverse 
temperature and L refers to the spatial box size. Further 
(Wf) with i G {1, 2} is the winding number squared in i- 
direction. By carefully investigating the spatial volumes 
and the J'/J dependence of p S iL, one can determine the 



III. DETERMINATION OF THE CRITICAL 
POINT AND THE CRITICAL EXPONENT v 

To calculate the relevant critical exponent v and to 
determine the location of the critical point in the pa- 
rameter space J' / J, one useful technique is to study the 
finite-size scaling of certain observablcs. For example, if 
the transition is second order, then near the transition, 
the observable p S iL p for i 6 {1,2} should be described 
well by the following finite-size scaling ansatz 

LP (t) = (l-b(LP)-n9o(t(LP) 1/ n, (3) 

where Olp stands for p S iL p , L p is the physical linear 
length of the system, t = {j c —j)/jc with j = (J'/J), b is 
some constant, v is the critical exponent corresponding 
to the correlation length £ and w is the confluent correc- 
tion exponent. Finally go appearing above is a smooth 
function of the variable ^ZJ 1 ) 1 /^. In practice, the LP ap- 
pearing in Eq. (J3J is conventionally replaced by the box 
size L used in the simulations when performing finite-size 
scaling analysis. We will adopt this conventional strategy 
in the first part of our analysis as well. From Eq. ([3]) , one 
concludes that the curves of different L for Ol, as func- 
tions of J'/J, should have the tendency to intersect at 
critical point (J'/ J) c for large L. To calculate the critical 
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FIG. 3: Fits of p E iL (upper panel) and p a L (lower panel) 
to Eq. (J3J). While the circles and squares on these two panels 
are the numerical Monte Carlo data from the simulations, the 
solid curves are obtained by using the results from the fits. 



exponent v and the critical point (J /J) c , in the follow- 
ing we will apply the finite-size scaling formula, Eq. ([3]), 
to both p s \L and p S 2L. Without losing the generality, in 
our simulations we have fixed J to be 1.0 and have varied 
J'. Further, the box size used in the simulations ranges 
from L = 6 to L = 64. We also use large enough /3 so that 
the observables studied here take their zero-temperature 
values. Figure [2] shows the Monte Carlo data of p s L 
and p 8 iL as functions of the parameter J'/J. The figure 
clearly indicates the phase transition is likely second or- 
der since different L curves for both p s L and p S 2L tend 
to intersect at a particular point in the parameter space 
J'/J. What is the most striking observation from our 
results is that the observable p s L receives a much severe 
correction than p s2 L does. This can be understood from 
the trend of the crossing among these curves of different 
L in figure [H Therefore one expects a better determina- 
tion of v can be obtained by applying finite-size scaling 
analysis to p S 2L. Before presenting our results, we would 
like to point out that since data from large volumes might 
be essential in order to determine the critical exponent v 
accurately as suggested in [l]|, we will use the strategy 
employed in [l3| for our data analysis as well. A Taylor 
expansion of Eq. ([3]) up to fourth order in tL 1 ^ is used to 
fit the data of p S 2L. The critical exponent v and critical 



point (J'/J) c calculated from the fit using all the avail- 
able data of p s2 L arc given by 0.6934(13) and 2.51962(4), 
respectively. The upper panel of figure |3j demonstrates 
the result of the fit. Notice both v and ( J' / J) c we obtain 
are consistent with the corresponding results found in [l[ . 
By eliminating some data points of small L, we can reach 
a value of 0.700(3) for v by fitting p S 2L with L > 26 to 
Eq. ([3]). On the other hand, with the same range of L 
(L > 26), a fit of p s L to Eq. © leads to v = 0.688(2) 
and (J'/J) c = 2.5193(2), both of which are consistent 
with those obtained in [1| as well (lower panel in figure 
■ By eliminating more data points of p s L with small 
L, the values for v and (J'/J) c calculated from the fits 
are always consistent with those quoted above. What 
we have shown clearly indicates that one would suffer 
the least correction by considering the finite-size scaling 
of the observable p S 2L. As a result, it is likely one can 
reach a value for v consistent with the 0(3) prediction, 
namely v = 0.7112(5) if large volume data points for 
p S 2 are available. Here we do not attempt to carry out 
such task of obtaining data for L > 64. Instead, we 
employ the technique of fixing the ratio of spatial wind- 
ing numbers squared in the simulations. Surprisingly, 
combining this new idea and finite-size scaling analysis, 
even from the observable p s \L which is found to receive 
the most severe correction among the observables con- 
sidered here, we reach a value for the critical exponent 
v consistent with v = 0.7112(5) without additionally ob- 
taining data points for L > 64. The motivation behind 
the idea of fixing the ratio of spatial winding numbers 
squared in the simulations is as follows. First of all, as 
we already mentioned earlier that the box size L used 
in the simulations is conventionally used as the L p in 
Eq. (J3j) when carrying out finite-size scaling analysis. 
For isotropic systems, such strategy is no problem. How- 
ever, for anisotropic cases, the validity of this common 
wisdom of treating L as LP is not clear. In particular, 
the same L does not stand for the same LP of the sys- 
tem for two different anisotropics J'/J. Hence one needs 
to find a physical quantity which can really character- 
ize the physical linear length of the system. Secondly, 
in magnon chiral perturbation theory which is the low- 
energy effective field theory for spin- 1/2 antiferromagnets 
with O(N) symmetry, an exactly cubical space-time box 
is met when the condition /3c = L is satisfied, here c is 
the spin- wave velocity and /3, L are the inverse temper- 
ature and box size as before. For spin-1/2 XY model 
on the square lattice, for large box size L, the numerical 
value of c determined by L//3 using the /3 with which 
one obtains (W 2 ) = 1/2({W?) + (IF 2 2 )) = (VF t 2 ) in the 
Monte Carlo simulations agrees quantitatively with the 
known results in the literature [14|. This result implies 
that the squares of winding numbers are more physical 
than the box sizes since an exactly cubical space-time 
box is reached when the squares of spatial and tempo- 
ral winding numbers are tuned to be the same in the 
Monte Carlo simulations. Consequently the physical lin- 
ear lengths of the system should be characterized by the 
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squares of winding numbers, not the box sizes used in 
the simulations. Based on what we have argued, it is 
(W?)/{Wi), not (L 2 /Li) 2 , plays the role of the quan- 
tity {L2/ L\) 2 for the system, here again we refer L\ with 
i G {1,2} as the physical linear length of the system in i- 
direction. As a result, fixing the ratio of spatial winding 
numbers squared in the simulations corresponds to the 
situation that the physical shape of the system remains 
fixed in all calculations. Indeed it is demonstrated in 
Q that rectangular lattice is more suitable than square 
lattice for studying the spatially anisotropic Heisenberg 
model with different antiferromagnetic couplings J\, J2 
in 1- and 2-directions. The idea of fixing the ratio of 
spatial winding numbers squared quantifies the method 
used in Q. 

The method of fixing the ratio of spatial winding num- 
bers squared is employed as follows. First of all, we per- 
form a trial simulation to determine a fixed value for the 
ratio of spatial winding numbers squared which we de- 
note by Wf and will be used later in all calculations. Sec- 
ondly, instead of fixing the aspect ratio of box sizes L\ 
and L2 in the simulations as in conventional finite-size 
scaling studies, we vary the variables L\, L2 and J'/J 
in order to satisfy the condition of a fixed ratio of spa- 
tial winding numbers squared. This step involves a con- 
trolled interpolation on the raw data points. In practice, 
for a fixed L2 one performs simulations for a sequence 

Li = L2, L2 ± 2, L<z ± 4, The criterion of a fixed 

ratio of spatial winding numbers squared is reached by 
tuning the parameter J2 / J\ and then carrying out a lin- 
ear interpolation based on (w/wf)^ 1 ^ 2 ^ for the desired 
observables, here w refers to the ratio of spatial wind- 
ing numbers squared of the data points other than the 
trial one. Notice since only the ratio of the physical lin- 
ear lengths squared is fixed, it is natural to use L2 in 
the finite-size scaling ansatz Eq. (|3]) both for the anal- 
ysis of p s \ and p S 2- The validity of this unconventional 
finite-size scaling method can be verified by considering 
the transition induced by dimerization for the Heisen- 
berg model with a ladder pattern anisotropic couplings 
(figureg]). For b ~ 0.95(22) in Eq. ©, we obtain a good 
data collapse for the observable (p s i)i n L p 1 (= (p s i)i n L2)- 
Above the subscript "in" means the data points are the 
interpolated one. To make sure that the step of inter- 
polation leads to accurate results, we have carried out 
several trial simulations and have confirmed that the in- 
terpolated data points are reliable as long as the ratio is 
kept small (table 1). On the other hand, for b ~ 1.30(18) 
in Eq. a good data collapse is also obtained for the 
observable p s \Li, here p s \ arc the raw data determined 
from the simulations directly. Figure [5] shows a compari- 
son between the data collapse obtained by using the new 
unconventional method introduced above (upper panel) 
and by the conventional method (lower panel). For ob- 
taining figureEl we have fixed v = 0.7112, ui = 0.78, and 
(J/J') c = 0.52367, which are the established values for 
these quantities. As one sees in figure [SJ the quality of 
the data collapse obtained with the new method is bet- 



j'/J Ll L2 Wf/w (psl)in £sl 



0.53 96 96 0.9558(33) 


0.008188(22) 


0.008198(7) 


0.53 96 94 0.9549(32) 


0.008391(21) 


0.0084098(74) 


0.545 90 94 0.9594(35) 


0.016862(33) 


0.016835(15) 


0.545 90 90 0.9539(36) 


0.017651(35) 


0.017676(17) 


0.535 98 98 0.9591(28) 


0.011707(28) 


0.0117297(124) 


0.54 96 96 0.9592(29) 


0.014838(37) 


0.014846(13) 


0.525 96 96 0.9503(41) 0.0072255(225) 


0.0072579(66) 



TABLE I: Comparison between interpolated and original val- 
ues of p„i for several data points. The data points which are 
used for interpolation are obtained from the simulations with 
Li x (1/2 + 2) (except the last row which is obtained from a 
simulation with (Li +2) x L^)- The inverse temperature /3 
for these data points are fixed to /3J = 800. 
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FIG. 4: Heisenberg model with a ladder pattern of anisotropy. 



ter than the one obtained with the conventional method, 
thus confirming the validity of the idea to fix the ratio of 
winding numbers squared in order to studying the critical 
theory of a second order phase transition. 

As demonstrated above, in general for a fixed L2, one 
can vary L\ and J'/J in order to reach the criterion of 
a fixed aspect-ratio of spatial winding numbers squared 
in the simulations. For our study here, without obtain- 
ing additional data, we proceed as follows. First of all, 
we calculate the ratio (W 2 ) / (IF|) f° r the data point at 
J'/J = 2.5196 and L — 40 which we denote by Wf. We 
further choose L\ = L in our data analysis. After obtain- 
ing this number, a linear interpolation for p s \ of other 
data points based on (w/wf)^ 1 / 2 * 1 is performed in order 
to reach the criterion of a fixed ratio of spatial winding 
numbers squared in the simulations. The w appearing 
above is again the corresponding (W 2 )/(W2) of other 
data points. Here a controlled interpolation similar to 
what we have done in studying the ladder anisotropic 
Heisenberg model is performed as well. Further, since 
large volumes data is essential for a quick convergence of 
v as suggested in we make sure the set of interpo- 
lated data chosen for finite-size scaling analysis contains 
sufficiently many points from large volumes as long as 
the interpolated results are reliable. A fit of the inter- 
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FIG. 5: Comparison between the results of data collapse us- 
ing the new unconventional finite-size scaling method (up- 
per panel) described in the text and the conventional method 
(lower panel) for the ladder anisotropic Heisenberg model. 
The result in the upper panel is obtained from simulations 
with box sizes [L-x — 6) x Li, (L2 — 4) x L2, (L2 + 4) x £2 
for various values of J2/J1 if the interpolations from such 
simulations are reliable. 



polated (psi)inL data to Eq. ([3]) with oj being fixed to 
its 0(3) value (w = 0.78) leads to v = 0.706(7) and 
{J/J) c = 2.5196(1) for 36 < L < 64 (figure©. Letting to 
be a fit parameter results in consistent v = 0.707(8) and 
(J' I J) c = 2.5196(7). Further, we always arrive at consis- 
tent results with v = 0.706(7) and (J'/J) c = 2.5196(1) 
from the fits using L > 36 data. The value of v we 
calculate from the fit is in good agreement with the ex- 
pected 0(3) value v = 0.7112(5). The critical point 
(J' / J)c = 2.5196(1) is consistent with that found in [l[ 
as well. To avoid any bias, we perform another analysis 
for the raw p s \L data with the same range of L and J' / J 
as we did for the interpolated data. By fitting this set 
of original data points to Eq. (j3]) with a fixed lu = 0.78, 
we arrive at v = 0.688(7) and (J'/J) c = 2.5197(1) (fig- 
ure [6]), both of which again agree quantitatively with 
those determined in [![. Similarly, applying this uncon- 
ventional finite-size scaling to p S 2 would lead to a nu- 
merical value of v consistent with v = 0.7112(5). For 
instance, the v determined by fitting (p S 2)inL to Eq. (J3]) 
is found to be v = 0.706(7), which agrees quantitatively 



with the predicted 0(3) value (figure [8]). Finally we 
would like to make a comment regarding the choice of 
Wf. In principle one can use Wf determined from any 
L and from any J'/J close to (J'/J) c . However it will 
be desirable to choose Wf such that the set of interpo- 
lated data used for analysis includes as many data points 
from large volumes as possible. Using the w / obtained at 
J'/J = 2.5191 (J'/J = 2.5196) with L = 40 (L = 44), we 
reach the results of v = 0.704(7) and (J'/J) c = 2.5196(1) 
(v = 0.705(7) and {J'/J) c = 2.5196(1)) from the fit with 
a fixed cj = 0.78. These values for v and ( J' / J) c agree 
with what we have obtained earlier. Indeed as we will 
demonstrate in another investigation, the critical expo- 
nent v determined by the idea of fixing the ratio of spatial 
winding number squared in the simulations is indepen- 
dence of the chosen reference point. 



IV. DISCUSSION AND CONCLUSION 

In this paper, we revisit the phase transition driven 
by dimerization for the spin- 1/2 Heisenberg model with 
a spatially staggered anisotropy on the square lattice. 
We find that the observable p S 2L suffers a much less 
severe correction compared to that of p s iL, hence is 
a better quantity for finite-size scaling analysis. Fur- 
ther, we propose an unconventional finite-size scaling 
method, namely we fix the ratio of spatial winding num- 
bers squared. As a result, the physical shape of the sys- 
tem remains fixed in all simulations and analysis. With 
this new strategy, we arrive at v = 0.706(7) for the criti- 
cal exponent v which is consistent with the most accurate 
Monte Carlo 0(3) result v = 0.7112(5) by using only up 
to L = 64 data points derived from both p s \L and p S 2L- 
Interestingly, the x 2 /d.o.f. obtained from the fits using 
the interpolated data are better than those resulted from 
the fits using the raw data (figures [U [7] and [8]) . This ob- 
servation provides another evidence to support the quan- 
titative correctness of the new unconventional finite-size 
scaling we proposed here. 

It seems that when carrying out the finite-size scaling 
analysis for the observables considered here, the use of 
physical linear lengths of the system, which are charater- 
ized by the spatial winding numbers squared, would lead 
to a faster convergence of v. It will be interesting to apply 
a similar technique to other observables such as Binder 
cumulants as well. However, for Binder cumulants, the 
correction from interpolation will cancel out because of 
the definition of these observables. Therefore to further 
test the philosophy behind the unconventional finite-size 
scaling method proposed here might require some new 
ideas. Nevertheless, with this new unconventional finite- 
size scaling method, we have successfully solved the puz- 
zle raised in [l[ by showing that the anisotropy driven 
phase transition for the spin- 1/2 Heisenberg model with 
a staggered spatial anisotropy indeed belongs to the 0(3) 
universality class. Of course, the conventional finite-size 
scaling analysis is more convenient since one does not 
need to carry out interpolation on the raw data. How- 
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FIG. 6: Fit of interpolated (p s i) in L data to Eq. ©. While 
the circles are the numerical Monte Carlo data from the sim- 
ulations, the solid curves are obtained by using the results 
from the fit. 



ever for the subtle phase transition considered in this 
study, without obtaining data of gigantic lattices, a new 
idea which is more physical oriented such as the one pre- 
sented here is necessary. Still, to clarify the puzzle of an 
unconventional phase transition for the model studied 
here as observed in [l[ by simulating larger lattices and 
using the conventional finite-size scaling method is desir- 
able. However, such investigation is beyond the scope of 
our study. 
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FIG. 7: Fit of original p s \L data to Eq. (J3j) - While the circles 
are the numerical Monte Carlo data from the simulations, the 
solid curves are obtained by using the results from the fit. 
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FIG. 8: Fit of interpolated (p s2 )i n i data to Eq. ©. While 
the circles are the numerical Monte Carlo data from the sim- 
ulations, the solid curves are obtained by using the results 
from the fit. 



